Eigenvalue problems, spectral parameter power 
series, and modern applications 



Kira V. Khmelnytskaya^, Vladislav V. Kravchenko^, 

Haret C. Rosu^ 

"'"Faculty of Engineering, Autonomous University of Queretaro, Queretaro, Mexico 
^Department of Mathematics, CINVESTAV del IPN, Campus Queretaro, 
Apartado Postal 1-798, Arteaga # 5, Col. Centro, Queretaro, Qro. 

76001 Mexico, kravclienko@qro.cinvestav.mx 
IPICyT, Instituto Potosino de Investigacion Cientifica y Tecnologica, 
Apdo Postal 3-74 Tangamanga, 78231 San Luis Potosi, Mexico, hcr@ipicyt.edu. mx 



Abstract 

Our revie"w is dedicated to a wide class of spectral and transmission problems 
arising in different branches of applied physics. One of the main difficulties in 
studying and solving eigenvalue problems for operators with variable coefficients 
consists in obtaining a corresponding dispersion relation or characteristic equa- 
tion of the problem in a sufficiently explicit form. Solutions of the dispersion 
relation are the eigenvalues of the problem. When the dispersion relation is 
known the eigenvalues are found numerically even for relatively simple problems 
with constant coefficients because even in those cases as a rule the dispersion 
relation represents a transcendental equation the exact solutions of which are 
unknown. 

In the present review we deal with the recently introduced method of spec- 
tral parameter power series (SPPS) and show how its application leads to an 
explicit form of the characteristic equation for different eigenvalue problems in- 
volving Sturm-Liouville equations with variable coefficients. We consider Sturm- 
Liouville problems on finite intervals; problems with periodic potentials involving 
the construction of Hill's discriminant and Floquet-Bloch solutions; quantum- 
mechanical spectral and transmission problems as well as the eigenvalue prob- 
lems for the Zakharov-Shabat system. In all these cases we obtain a charac- 
teristic equation of the problem which in fact reduces to finding zeros of an 
analytic function given by its Taylor series. We illustrate the application of the 
method with several numerical examples which show that at present the SPPS 
method is the easiest in the implementation, the most accurate and efficient. 
We emphasize that the SPPS method is not a purely numerical technique. It 
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gives an analytical representation both for the solution and for the characteristic 
equation of the problem. This representation can be approximated by different 
numerical techniques and for practical purposes constitutes a powerful numer- 
ical method but most importantly it offers additional insight into the spectral 
and transmission problems. 

Keywords: spectral parameter power series; Sturm-Liouville problems; dispersion rela- 
tions; periodic potentials; Hill's discriminant; supersymmetry; Zakharov-Shabat sys- 
tem 

1 Introduction 

Solution of second-order linear differential equations belongs to a classical field of 
mathematics in which as an overwhelming majority of its users from students to active 
researchers believe that everything that was possible to do theoretically is essentially 
done, and whatever the analysts invent, in the domain of practical solution of equations 
and related models anyway the numerical discrete schemes more and more refined due 
to the massive efforts of experts in numerical analysis and computer sciences will be 
more accurate and efficient. For the good luck of mathematics and its applications a 
considerable number of analysts as well as of advanced users of mathematical methods 
in physics are aware of many strong limitations in applicability of available numerical 
schemes. For examples, if the discrete spectrum of a problem is not necessarily real, 
practically the whole machinery of advanced numerical techniques does not apply 
bringing to the surface in fact one option only: finite differences. The importance 
of this technique lies in its universality. Nevertheless it usually gives way to many 
other approaches whenever they become applicable. The main difficulty in finding 
complex eigenvalues is due to the fact that the universally used shooting method 
works if only there is a clear criterion for choosing every next shot. Meanwhile the 
real numbers is an ordered set and zero is located always between a negative and a 
positive outcomes of the corresponding shots, the complex plane does not admit such 
a simple rule. There are many other situations (even when the eigenvalues are real) 
when the shooting procedure finds considerable difficulties and at the same time the 
method of finite differences is not applicable at all. For example, when the spectral 
parameter participates in the boundary conditions. Such situation in fact is more 
common in applications than otherwise. 

In the present review we discuss an approach developed in the last few years and 
called the spectral parameter power series (SPPS) method. It is important to notice 
that the SPPS method is not merely another numerical technique. On the contrary, it 
is an analytical approach giving new analytical results and at the same time lending 
itself to numerical calculation. The SPPS method allows one to obtain two linearly 
independent solutions of the Sturm-Liouville equation (in Section |2] we specify the 
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conditions imposed on the coefficients) 

{pu'y + qu = Xru (1) 

in the form 

oo oo 

Ui{x) = ^^afc(x)A'^ and U2{x) = hk{x)X^ 

k=0 k=0 

where A is a spectral parameter and the series are uniformly convergent. Such repre- 
sentations from time to time appear in mathematical literature in different contexts. 
We mention here |6l Sect. 10] and [28j. The main difference is the form in which the 
coefficients and b^, k = 0,1, . . . are represented. In previous works the calculation 
of the coefficients was proposed in terms of successive integrals with the kernels in 
the form of iterated Green functions (see [HI Sect. 10]). This makes any computation 
based on such representation difficult and less practical. Moreover, theoretical study 
of the corresponding series and their properties becomes considerably more compli- 
cated. We show that a^, and bj. can be calculated in terms of the coefficients p and r 
of ([T]) as well as of a particular solution of the equation {pv')' + qv = 0. The obtained 
representations of the coefficients and bk are relatively simple and well suited both 
for theoretical estimates and for numerical computation with a minimum of program- 
mer's efforts required. Behind this representation of the coefficients and bk there is 
one of the possible factorizations of the operator L = -^p-^ + q sometimes called the 
Polya factorization (see [53]). 

The main advantage of the SPPS method is akin to that of such asymptotic ap- 
proaches as the WKB method - it allows one to work with an analytical representation 
of the solution instead of a table of values delivered by a numerical method. This is 
important from different points of view, often it gives a new physical insight into the 
problem. The difference between the SPPS and the asymptotic techniques lies in the 
fact that to apply the SPPS method it is not necessary to assume the smallness or the 
largeness of the parameter A. For example, the WKB approximation can be efficiently 
applied to ([1]) when A is sufficiently large which is quite useless when an eigenvalue 
problem related to ([1]) is considered. 

We show that 1) for solving initial value and boundary value problems the SPPS 
method performs better or equal in comparison to purely numerical techniques; 2) it 
is highly advantageous when the solution is required for many different values of the 
spectral parameter; 3) the SPPS method allows us to write down an explicit form of the 
characteristic equations for many different spectral problems which in practice reduces 
the spectral problem to finding zeros of a corresponding analytic function given by its 
Taylor series. We emphasize that the method is applicable in different situations when 
some other approaches are unavailable (complex eigenvalues; A-dependent boundary 
conditions, etc.) The method is simple and can be introduced in mathematical courses 
for physicists. 

In Section |2l we review the main results from [62] and [66] concerning the SPPS 
representation of the solutions of ([T]) and show that even in solving initial and boundary 
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value problems this technique converted into a simple numerical algorithm is clearly 
competitive when compared to standard routines for numerical integration of linear 
ordinary differential equations. In Section [3] we apply the SPPS method to Sturm- 
Liouville spectral problems with or without the spectral parameter in the boundary 
conditions. Here together with some results from [66] we present new results con- 
cerning the problems admitting complex eigenvalues. Section H] is dedicated to the 
spectral problems for periodic potentials. We give an SPPS representation for Hill's 
discriminant [55] and show how the SPPS method allows one to construct the Bloch 
solutions of the problem. In Section [5] we consider two classical problems of mathe- 
matical physics: the quantum-mechanical spectral problem and the transmission prob- 
lem. Following [IH] we present a characteristic (dispersion) equation equivalent to the 
eigenvalue problem for the Schrodinger operator with a potential which is an arbitrary 
continuous function on a finite interval outside of which it is constant. We discuss 
some numerical tests as well. The transmission problem is presented in the context of 
electromagnetic wave propagation as the problem of calculation of the reflection and 
transmission coefficients for a plane wave which is incident on an inhomogeneous layer 
under an arbitrary angle of incidence. Following [T7] we discuss application of the 
SPPS method to this problem. Section E] is dedicated to the eigenvalue problem for 
the Zakharov-Shabat system. We present the dispersion equation [68] for the prob- 
lem for a real-valued, finitely supported potential and discuss its practical application. 
Finally, in Section [7] we make some concluding remarks. 

This review is for the colleagues interested in all kinds of problems involving so- 
lution of Sturm-Liouville type equations. We hope to attract more attention to the 
SPPS approach which combines the possibility to work with analytical representa- 
tions of solutions and of characteristic equations of the problems with simplicity, rapid 
convergence, and accuracy when used for numerical computation. 

2 Spectral parameter power series representation 
for solutions of the Sturm-Liouville equation 

Let us consider the Sturm-Liouville equation 



where p, q and r are complex valued functions and A is a complex parameter. The 
following result [156] gives us a convenient form for a general solution of ([2]) as a spectral 
parameter power series. 

Theorem 1 Assume that on a finite interval [a, 6], equation 




(2) 




(3) 
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possesses a particular solution Uq such that the functions UqT and 1/{uqp) are contin- 
uous on [a,b]. Then the general solution of ([^ on {a,b) has the form 



U = CiUi + C2U2 , 

where c\ and 02 are arbitrary complex constants, 

00 00 
u^ = uoJ2^'x'^''^ and U2 = uoJ2^'X^''+''> 

k=0 k=0 

with X*^"^ and X^"'^ being defined by the recursive relations 

^(0) ^ 1 ^(0) ^ T 



(4) 



(5) 



/ X 



X) 



X("-i)(s)Mg(s)r(s)rfs, n odd, 



^Hs)^7^r-r^ ds, n even, 



(6) 



(7) 



( X 



-2,] , , ds, n odd. 



X) = < 



X*^" ^\s)u1{s)r{s) ds, n even. 



XO 



where Xq is an arbitrary point in [a,b] such that p is continuous at Xq and p{xq) 7^ 0. 
Further, both series in ^ converge uniformly on [a,h\. 

For a detailed proof we refer to [66j . It is based on some simple observations. First 
of all the knowledge of a particular solution of ([3]) allows one to factorize the Sturm- 
Liouville operator L = -^p-j^+q in the form L = vui-^ — , also known as the Polya 

^ ax^ ax ^ MO dx ^ ^>ax uq' 

factorization [53] as mentioned before. This form is well suited for establishing how 
the operator acts on each member of the series ([5]). For example, ^uqX^^'^^^ = 
^^^{2fc-2)^ A; G N. Analogously, \L {uoX^'^''+^^) = mqX^^'^-i). 

Considering the system of functions {ipn}'^=o defined as follows ipo = uq, 



(9) 



Mo(x)X(")(x), n odd, 
Uo{x)X^"\x), n even, 

we find that ^Lipo^i = and ^Lipn = n = 2,3, . . .. These properties are char- 

acteristic for the so-called L-bases (introduced in [M], unfortunately this important 
book has not been translated into English), and hence formulas ([S])-® represent a 
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practical way to calculate an L-basis. In [M] it was established that the system of 
functions {fn}'^=Q is complete in ^2(0, 6). For further related properties we refer to 

m- 

To establish the uniform convergence of the series in ([5]) as well as to get a rough 
but useful estimate for the velocity of their convergence it is sufficient to observe that 



X^^"^ < {max\rul\f( 



max 



fc 1, |2fc 

— a 



(10) 



00 



V puo J i2k)\ 

(a similar inequality is available for |X*^^'^'+^'' | as well). Thus, the series IS 

fc=0 

00 

majorized by a convergent numerical series (^yy with c = |A| (max \ruQ\) ^max 

k=0 

It is worth noticing that from ([5]) it is easy to obtain that ui and U2 satisfy the 
following initial conditions 

Ml(Xo) = Uo{Xo), u[{Xo) = Uq{Xo), (11) 

U2{xo) = 0, u'^ixo) = f \ f (12) 

Remark 2 The possibility mentioned in the Introduction to represent solutions of the 
Sturm- Liouville equation in the form of spectral parameter power series is by no means 
a novelty, though it is not a widely used tool. In fact, besides the work reviewed below 
we are able to mention only Sect. 10], J2^ and the recent paper J6J^) and to the 
best of our knowledge it was applied for the first time for solving spectral problems in 
l6^ . The reason of this underuse of the SPPS lies in the form in which the expansion 
coefficients were sought. Indeed, in previous works the calculation of coefficients was 
proposed in terms of successive integrals with the kernels in the form of iterated Green 
functions (see Sect. 10]). However, this makes any computation based on such 
representation difficult, less practical, and even proofs of the most basic results like, 
e.g., the uniform convergence of the spectral parameter power series for any value of 
A e C (established in TheoremU\) are not an easy task. For example, in p. 16] the 
parameter A is assumed to be small and no proof of convergence is given. Moreover, in 
\TT] a discrete analogue of TheoremUltogether with some further applications to Jacobi 
operators were established and it was pointed out that as well as in the continuous case 
the SPPS representation for solutions of the Jacobi operators was considered as a 
perturbation technique, however, even the situation with the convergence of such series 
was not satisfactorily understood. We recommend the book ^ where the possibility of 
divergence of such series as those considered in the present work is assumed. Due to 
the representation of the expansion coefficients similar to it was shown in [77]/ 

that the series are not only convergent but in the discrete case they are actually finite 
sums. 



6 



The way of how the expansion coefficients in ^ are calculated according to Ij^ 
and ^ is relatively simple and straightforward. This is why the estimation of the 
rate of convergence of the series ^ presents no difficulty, see / fiOj) . Another crucial 
feature of the introduced representation of the expansion coefficients in ^ consists in 
the fact that as we repeatedly observe in subsequent pages not only the expansion coef- 
ficients themselves also denoted by in ^ are necessary in solving different spectral 
problems related to the Sturm- Liouville equation but also the formal powers apparently 
obtained as a by-product of the recursive integration procedure ^ and namely the 
functions X*^^^"*"^) and X^'^^\ A; = 0, 1, 2, . . which do not participate explicitly in the 
representation of solutions naturally appear in dispersion equations corresponding 
to the spectral problems. See / f^) and l[2^) for the dispersion equations equivalent 
to the classical Sturm- Liouville eigenvalue problem, / f5^) for the SPPS representation 
of the Hill discriminant, / f55]) for the dispersion equation equivalent to the quantum- 
mechanical eigenvalue problem on the whole axis or ( fg3]) for the dispersion equation 
equivalent to the Zakharov-Shabat eigenvalue problem. 

The consideration of formal powers ^ and ^ as infinite families of functions 
intimately related to the corresponding Sturm- Liouville operator led in fW^ . fl^ and 
^U7j to a deeper understanding of the transmutation operators f^, [T^ also known 
as transformation operators f7g| /, J77| /. Indeed, the functions (pn{x) resulted to be the 
images of the powers under the action of a corresponding transmutation operator 
fT3f . This makes it possible to apply the transmutation operator even when the oper- 
ator itself is unknown (and this is the usual situation - there are very few explicitly 
constructed examples available) due to the fact that its action on any polynomial is 
known. This result was used in IW^ and fT3f to prove the completeness (Runge-type 
approximation theorems) for families of solutions of two-dimensional Schrddinger and 
Dirac equations with variable complex-valued coefficients. 

Remark 3 One of the functions ru^ or 1/ {pu^ may not be continuous on [a, h] and yet 
Ui oru2 may make sense. For example, in the case of the Bessel equation {xu')' — -u = 
—Xxu, we can choose uo{x) = x/2. Then l/{pul) ^ C[0, 1]. Nevertheless, all integrals 

Ij^ exist and ui coincides with the nonsingular Ji(a/Ax), while U2 is a singular 

solution of the Bessel equation. 

Remark 4 When p and q are real-valued, p{x) ^ for all x G [a, h] and p, p' , q are 
continuous functions on [a,b], the equation 

Lv = (13) 

is a regular Sturm- Liouville equation and possesses two linearly independent real valued 
solutions Vi and V2. Due to Sturm's separation theorem (see, e.g., /57, p. 10]) their 
zeros occur alternately and hence Uq = Vi-\-iv2 can be chosen as the required in theorem 
Ulparticular solution. Ifr is a continuous on [a,b] (in general, complex-valued) function 
then the conditions of theorem [I] are fulfilled. 
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The solutions Vi and V2 from Remark 4 can be in fact calculated using the same 
procedure from Theorem [TJ Indeed, consider equation fll3p which can be written in 
the form 

{pv'y = —qv. 

It has the form (|2]) with r := — g, A = 1 and with a convenient solution f = 1 of the 
(homogeneous) equation {pv'q)' = 0. Apphcation of theorem [T] gives us the following 
two hnearly independent solutions of (IT^ . 



where 



k=0 n=0 



yiO) = 1^ y(0) ^ 1^ 

nodd, 
n even, 



(14) 



(15) 



(16) 



( X 



F(")(a;) = < 



nodd. 



(17) 



even, 



and the series in the equalities for vi and V2 converge uniformly on [a, h]. 
Note that 



fi(xo) 



v[{xq) = 0, t;2(a;o) = 0, v'^ixo) = l/p{xo). 



Solutions of ( IT3I) in the form f lT^ is a long known result (see, e.g., |100] ). 

Before we proceed to discuss eigenvalue and scattering problems it is worth noticing 
that the representation of a general solution of the Sturm-Liouville equation in the form 
of a spectral parameter power series (SPPS) given by theorem [1] represents a natural 
and highly competitive method for numerical solution of initial and boundary value 
problems. Compared to the best standard routines it performs better or equal and 
with minimal programmer's efforts. Moreover, the advantages of using SPPS become 
even more transparent when the solution of the problem is required for many different 
values of the spectral parameter. In such case the auxiliary functions X^"'^ and X*^"\ 
n = 0, 1, 2, . . . should be computed only once and then substitution of values of A into 
the expressions gives us a solution of equation ([2]) for as many different values of the 
spectral parameter as needed at no additional computational cost. Nevertheless, first. 
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let us show how SPPS performs at the terrain of numerical ODE solvers for solution 
of initial value problems. In [17j we made use of Matlab 7 and as a first step compared 
our results with standard Matlab ODE solvers [3] , [92] , especially with ode45 which in 
the considered examples gave always better results than other similar programs. Here 
we give two examples from |T7]. 

For numerical approximations we consider partial sums of the infinite series (jS]) 

N N 

and ( Il4|) . e.g., ui = mq^A'^X^^'') and U2 = mo^A''X(^''+^). The algorithm was 

k=0 k=0 

implemented in MATLAB. For the recursive integration we have chosen the following 
strategy. On each step the integrand is represented through a cubic spline using the 
spapi routine and the integration is performed using the fnint routine (both from the 
spline toolbox of MATLAB). 

Consider the following initial value problem for f|T3|) : p = 1, q = c^, v{0) = 1, 
v'(0) = —1 on the interval (0, 1). For c = 1 the absolute error of the result calculated 
by ode45 (with an optimal tolerance chosen) was of order 10~^ and the relative error 
was of order 10~^ whereas the absolute error of the result calculated with the aid of 
the SPPS representation with N from 55 to 58 was of order 10~^^ and the relative 
error was of order 10"^'^. Taking c = 10 under the same conditions the absolute and 
the relative errors of ode45 were of order 10~^ and 10~^ respectively meanwhile our 
algorithm gave values of order 10~^^ in both cases. For the initial value problem: 
p = —1, q = (?, f (0) = 1, v'(0) = —1 on the interval (0, 1) in the case c = 1 the 
absolute and the relative errors of ode45 were of order 10~^ whereas in our method 
this value was of order 10~^^ already for N = 50. For c = 10 the absolute and the 
relative errors of ode45 were of order 10~^ and 10~^ respectively and in the case of our 
method these values were of order 10~^^ and 10~^^ for N = 50. 

Consider yet another example. Let p = —1, q{x) = c^x^ + c in f lT^ . In this case 
the general solution has the form 

v{x) = e'="'/2 ^ ^2 e-'^'d^j . 

Take the same initial conditions as before, v{0) = 1, v'{0) = —1. Then, while for c = 1 
the absolute and the relative error of ode45 were both of order 10~* and for c = 30 the 
absolute error was 0.28 and the relative error was of order 10~^, our algorithm [N = 58) 
gave the absolute and relative errors of order 10~^^ for c = 1 and the absolute and 
relative errors of order 10~^ and 10~^^ respectively for c = 30. All calculations were 
performed on a common PC with the aid of Matlab 7. 

The results of our numerical experiments show that in fact the SPPS represen- 
tations offer a powerful method for numerical solution of initial value and boundary 
value problems for linear ordinary differential second-order equations. The numerical 
calculation of the involved integrals does not represent any considerable difficulty and 
can be done with a remarkable accuracy. 
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Another observation which can be of use in many different situations is that if for 
some purposes a derivative of the solution of ([2]) is required there is no need to apply 
to the obtained solution an algorithm for numerical differentiation. Instead, it is easy 
to see that 



OO 

Ml = Ui 

k=l k=0 



Thus, the calculated auxiliary functions X^"^ and X^''^\ n = 0, 1, 2, . . . are used once 
more, this time for obtaining the derivative of the solution. 



3 Solution of Sturm-Liouville problems 

In this section we outline the main ideas behind the application of the SPPS method 
to the solution of Sturm-Liouville eigenvalue problems referring the interested reader 
to for additional details and numerical examples. The SPPS method allows one 
to reduce the Sturm-Liouville problem to the problem of finding zeros of an analytic 
function of the complex variable A. Numerically the problem is reduced to finding roots 
of a polynomial in A. To find the precise expressions for Taylor coefficients of that 
analytic function let us consider the general Sturm-Liouville problem with unmixed 
boundary conditions. Thus, we look for the eigenvalues and eigenfunctions of the 
problem 

{pu'y + qu = Xru, (19) 

'u(a) cos a + «'(«) sin a = 0, (20) 

M(6)cos/3 + n'(6)sin/3 = , (21) 

where [a, b] is a finite segment of the x-axis, a and f3 are arbitrary real numbers. 

Let us choose the point xq from theorem [1] being equal to a and consider the 
solutions ui and U2 of (|T9l) defined by ([5]). Then from (fTTl) and ( fT2l) we obtain that a 
linear combination u = ciUi + C2M2 satisfies the following conditions at a: 

u{a) = ciUo{a) and u{a) = ciUQ^a) + C2/{uo{a)p{a)). 

Thus, in order that u satisfy fl20l) . the constants Ci and C2 must satisfy the equation 

sin (y 

Ci(uo(a) cos a + u'n(a) sin a) + C2 — i—r~rT = O5 

Uo{a}p{a) 

which gives C2 = 7C1 when a 7^ nn, with 7 = —uo{a)p{a){uo{a) cot a + 'Ug(a)), whereas 
Ci = when a = Tcn. In the latter case we have that if an eigenfunction of the problem 
for a given A exists, up to a multiplicative constant it must have a form u = U2- The 
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second boundary condition f l2ip together with fll8l) leads to the following characteristic 
equation for the eigenvalues 

oo / oo ^ oo \ 

cos/3MoW$^A'=X(2'=+i)(6)+sin/3 u'Q{b)Y^X^X^^''+'\b) + | A ^A^X(^^)(&) = 

fc=0 V fc=0 ' k=0 J 

which is the same to 

T^' (x^''+'\b) (cos/3«o(6) + sin/3 w'o(&)) + T/^,, X(''\b)] = 0. (22) 
t^o \ uo{h)p(h) J 

Thus, the Sturm-Liouville problem (fT9|) -( l2T]) in the case a = nn reduces to find zeros 

oo 

of the analytic function ^^OfcA'^ where the Taylor coefficients ak have the form ak = 

(cos/3no(6) + sm°/3<(6)) + ^^X^^'^b). 
Now let us suppose a ^ vrn. Then the boundary condition (12T|) implies that 

/ oo oo \ 

(uo(6) COS/3 + u'o(6) sin/3) ^A'=X(2*^)(6) + 7 J]a'=X(2'=+i)(6) 



vfc=0 A;=0 



+ ;74!^ I + tE^'^^^'H^) ) = 0. (23) 



Mo(&)p(&) 

problem 

of the analytic function k(A) = ^^amA*" where 

m=0 

ao = (mo(6) cos/3 + m'o(6) sin/3) (1 + 7^^'H&)) + 



,fe=l fc=0 

Thus the spectral problem ([2]), (1201) . (pTj) reduces to the problem of calculating zeros 



Mo(6)p(6) 

and ^ 

a„ = (mo(6) cos/3 + m'o(6) sin/3) (x^''"'-\b) + 7X(2'"+i)(6) 



+ ^!^ (x(2-^)(6) + 7X(2-)(6)) , m = 1, 2, . . . . 

This reduction of a Sturm-Liouville spectral problem to finding zeros of an analytic 
function given by its Taylor series lends itself to a simple numerical implementation. To 

N 

calculate the first n eigenvalues we consider the Taylor polynomial KAf(A) = ^^amA"^ 

m=0 

with N > n. Thus the numerical approximation of eigenvalues of the Sturm-Liouville 
problem reduces to the calculation of roots of the polynomial /tAr(A). 

In many physical applications (see O |20l |23l |26l HH [99] and references therein) 
the Sturm-Liouville problems with boundary conditions dependent on the spectral 
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parameter arise. In this case together with equation (fT^ and boundary condition ( 12U]) 
the eigenfunction must satisfy a second boundary condition of the form 

- P2u{h) = ^{\) {P[u{h) - , (24) 

where is a complex- valued function of the variable A and /32, P'l, P'2 are complex 
numbers. For some special forms of the function tp such as (^(A) = A or 99(A) = A^ + 
C1A-I-C2, results were obtained [23], [12] concerning the regularity of the problem (IT^ . 
f l20|) . f l24|) : we will not dwell upon the details. In general, the presence of the spectral 
parameter in boundary conditions introduces additional considerable difficulties both 
in theoretical and numerical analysis of the problems. Nevertheless the SPPS approach 
gives a simple and natural insight into the problem, and its practical application for 
numerical calculations is not in fact more difficult than in the previously considered 
situation of A-independent boundary conditions. 

For simplicity, let us suppose that a = and hence the condition (120|) becomes 
u{a) = 0. Then as was shown above, if an eigenfunction exists it necessarily coincides 
with U2 up to a multiplicative constant. In this case condition becomes equivalent 
to the equality [5S] 

00 / \ 00 

{uo{b)^,{X) - <(6)^2(A)) J2x'^X(''^'\b) - = ' (25) 

where v5i,2(A) = /3i,2 — f^i 2fW- Calculation of eigenvalues given by ( l25l) is especially 
simple in the case of (p being a polynomial of A. Precisely this particular situation 
was considered in all of the above-mentioned references concerning Sturm-Liouville 
problems with spectral parameter dependent boundary conditions. In any case the 
knowledge of an explicit characteristic equation ( !25ll for the spectral problem (fT9l) . 
(1201) . (12^ makes possible its accurate and efficient solution. 

The paper [66] contains several numerical tests corresponding to a variety of com- 
putationally difficult problems. All they reveal an excellent performance of the SPPS 
method. We do not review them here referring the interested reader to [66]. Instead we 
consider another interesting example, a Sturm-Liouville problem admitting complex 
eigenvalues. 

Consider the equation ([2]) with p = —1, g = and r = 1 on the interval (0, tt) 
with the boundary conditions u{0) = and ^(Tr) = —X'^u{tt). The exact eigenvalues 
of the problem are A„ = together with the purely imaginary numbers X± = ±i. 
Application of the SPPS method with = 100 and 3000 interpolating points (used 
for representing the integrands as splines) delivered the following results Ai = 1, 
A2 = 4.0000000000007, A3 = 9.00000000001, A4 = 15.99999999996, A5 = 25.000000002, 
Ae = 35.99999997, A7 = 49.0000004, Ag = 63.9999994, Ag = 80.9996, Aio = 100.02 and 
X± = ±i. Thus, the complex eigenvalues are as easily and accurately detected by the 
SPPS method as the real eigenvalues. Note that for a better accuracy in calculation of 
higher eigenvalues of a Sturm-Liouville problem an additional simple shifting procedure 
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(described in [BS]) based on the representation of solutions of not as a series in 
powers of A but in powers of (A — Aq) is helpful. We did not apply it here and hence 
the accuracy of the calculated value of Aio is considerably worse than the accuracy 
of the first calculated eigenvalues which in general can be improved by means of the 
mentioned shifting procedure. 

4 Periodic potentials: Floquet-Bloch solutions and 
Hill's discriminant 

Towards the end of the 19th century, Hill |16] and Floquet [38] initiated the rigor- 
ous study of the spectral properties of periodic Sturm-Liouville equations with real 
coefficients and with periodic (and antiperiodic) conditions imposed to their solutions 
y(0) = ±y(7r), where y = {y,y')'^^. In 1883, Floquet established Floquet's theorem 
asserting that every solution is a linear superposition of independent solutions, both 
of the form of exponential factors multiplied by periodic functions, while Hill's work 
published in 1877 in Cambridge, Mass., and reprinted in Europe in 1886, deals with 
a special component of the motion of the lunar perigee. Other fundamental results 
concerning the sequence of eigenvalues have been obtained by Lyapunov in 1902 [73j . 

The paper of Hill made this class of equations of interest for many authors and the 
term Hill equations has been commonly used since about a century for second-order 
linear differential equations with periodic coefficients, in particular for y" + f{t)y = 
with f{t) a periodic function. There are three basic methods to determine the existence 
of eigenvalues in this case, (i) through Floquet theory and Hill's discriminant -D(A), (ii) 
using operator theory and variational techniques, and (iii) via Priifer's transformation. 
The standard method is the first one which essentially reduces itself to seeking the 
nontrivial solutions (also known as quasi-periodic solutions) of a Hill-type equation 
that satisfy 

y(vr) = /3y(0) 

for some complex parameters /3 known as 'Floquet multipliers' which solve the quadratic 
equation 

(3^ - D{X)I3 + 1 = . 

Thus, to get the Floquet multipliers one needs Hill's discriminant, which in general 
is a real- valued function expressed in terms of fundamental solutions of the given Hill 
equation. Notice that from the standpoint of the Floquet multipliers the periodic and 
antiperiodic solutions are special cases corresponding to (3p = e^*^'^* and = e^'^''^^^'^\ 
respectively, for any integer k. 

In the realm of quantum mechanics, the first fundamental application of Floquet 
theory belongs to Bloch |8] in 1928, who, using the quantum physical terminology, 
introduced Floquet's theorem for the special case in which the exponential factors are 
plane waves in the important context of solid state physics. Interestingly, no reference 
to the mathematical literature is given in Bloch's paper. In 1931, Kronig and Penney 
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[nS] introduced the first model of bands (stability regions) and gaps (instability regions) 
for the motion of electrons in crystal lattices. The activity in this area has been retaken 
only after the second world war period |5T1 |T6l |9ll [91] and at the present time a vigorous 
progress takes place in extended research lines covering nanostructures and photonic 
crystals. 

In the rest of this section, we briefiy describe our recent result [H] of expressing 
the Hill discriminant in terms of SPPS for a Hill-Sturm-Liouville equation of the type: 

- (p(x)/'(x, A))' + g(x)/(x, A) = Xfix, A) . (26) 

We assume that p{x) > 0, p'{x) and q{x) are continuous bounded periodic functions 
of period T. 

We first recall some necessary definitions and basic properties from the Floquet 
(Bloch) theory. For more details see, e.g., [76l [33] . 

For each A there exists a fundamental system of solutions, i.e., two linearly inde- 
pendent solutions of ( 12^ fi{x, A) and /2(x, A) which satisfy the initial conditions 

/i(0,A) = l, /((0,A) = 0, /2(0,A) = 0, /^(0,A) = 1. (27) 

The Hill discriminant has the following expression in terms of the fundamental system 
of solutions 

D(A) = /i(T,A) + /^(T,A). 

Employing -D(A) one can easily describe the spectrum of the corresponding equation. 
Namely, the values of A for which |-D(A)| < 2 form the allowed bands or stability 
intervals meanwhile the values of A such that |-D(A)| > 2 belong to forbidden bands or 
instability intervals [76]. The band edges (values of A such that |-D(A)| =2) represent 
the discrete spectrum of the operator, i.e., they are the eigenvalues of the operator 
with periodic (-D(A) = 2) or antiperiodic (-D(A) = —2) boundary conditions. The 
eigenvalues A„, n = 0, 1,2, ... form an infinite sequence Aq < Ai ^ A2 < A3..., and an 
important property of the minimal eigenvalue Aq is the existence of a corresponding 
periodic nodeless solution /o(x, Aq) [76]. In general solutions of ( 126|) are not of course 
periodic, and one of the important tasks related to Sturm-Liouville equations with 
periodic coefficients is the construction of quasiperiodic solutions. In this paper, we 
use the matching procedure from [STJ for which the main ingredient is the pair of 
solutions /i(x. A) and f2{x, A) of satisfying conditions ([27|) . Namely, using /i(x. A) 
and f2{x, A) one obtains the quasiperiodic solutions f±{x + T) = (3±f±{x) as follows 

/±(x,A) = /3^F±(x-nT,A), { ""Jf o!±lI ±2, .^.^^ ' ^^^^ 

where F±{x, A) are the so-called self-matching solutions, which are the following linear 
combinations F±{x,X) = fi{x,X) + a±f2{x,X) with a± being roots of the algebraic 
equation f2{T,XW + - f2{TA))a - /{(T,A) = 0. The Bloch factors /3± 

are a measure of the rate of increase (or decrease) in magnitude of the self-matching 



14 



solutions F±{x,X) when one goes from the left end of the cell to the right end, i.e., 
I3±{\) = '^^(o'^j • The values of I3± are directly related to the Hill discriminant, /3±(A) = 

1{D{X)t\/D^{X) -4), and obviously at the band edges /?+ = /?_ = ±1 for D{X) = ±2, 
correspondingly. 

4.1 The SPPS series representation of Hill's discriminant 

The SPPS construction method of the solutions A) and /2(x, A) satisfying the ini- 
tial conditions fl27|) is based on the knowledge of one non- vanishing particular solution 
/o(x, Ao) bounded on [0, T] together with j^^^^^- In the case of Hill's equation the first 
eigenvalue Aq generates a nodeless periodic eigenfunction /o(x,Ao). In what follows, 
we initially suppose that the value of Aq is known. Note that, it can be obtained by 
different methods including the same SPPS method [66j as we explain in subsection 

KB 

Given Aq, we proceed in three steps in order to obtain the representation of Hill's 
discriminant: 

• the first one is the construction of a particular nodeless solution fo{x, Aq) 
which is periodic, i.e., /o(x + T, Aq) = fo{x, Aq), 

• the second one is the construction of the fundamental system of solutions 

A) and f2{x, A) for all values of the parameter A, 

• the final step is getting the representation of Hill's discriminant. 
We detail each of the steps in the following subsections. 

4.2 The nodeless periodic solution 

We want to obtain the nodeless periodic solution /o(x, Aq) for A = Aq of the equation 

- (p(x)(/o(x))')' + g(x)/o(x) = Ao/o(x) . (29) 

To achieve this goal, we first have to construct in SPPS form the fundamental system 
of solutions of fl29|) . These solutions are not necessarily periodic. However, one can 
follow the old procedure of James |5T] allowing to obtain from /o,i(a;, Aq) and fo,2{x, Aq) 
the Floquet type solutions which degenerate to a single periodic/antiperiodic solution 
/o(x, Aq) since Aq represents a band edge. 

The functions /o,i(a^, Aq) and /o,2(2^, Aq) can be calculated according to iteration 
formulas of the type (fT5|) - (fT7|) 

oo oo 

/o,i(x,Ao)= Yl ^o"^ /o,2(x,Ao)=p(0) J] (30) 

even n=0 odd n=l 
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where 



- 1, ^ 1, 



/;x("-^)(0(g(0 - Ao)rfe for an odd n 

lo ^(S""^^ (0 ^c?^ for an even n 

/J" X^""^^ (0 ^c?^ for an odd n 



lo ^o" ^^(Ol-ylO - ^o)d^ for an even n . 

The periodic nodeless solution of ( l29l) is constructed as a particular case of a quasi- 
periodic solution ( |28l) . essentially as a self- matching solution, i.e., 



/o(x, Ao) = /o,i(2; - nT, Ao) + «p/o,2(a; - r^T, Aq), (31) 
riT < X < {n + 1)T 



n = 0,l,2,... , 

since the Floquet phase multiplier is /5 = 1 in the periodic case and ap = ■^o,2(^'^^o)^^(),i(^.-^o) ^ 
see 



4.3 Fundamental system of solutions 

Once having the function fo{x, Aq), the solutions /i(x. A) and /2(x, A) for all values of 
the parameter A can be given using the SPPS method once again 

/i(x. A) = ^So(a;, A, Aq) + p(0)/^(0)/o(x)Si(x, A, Ao), 

(32) 

Mx,X) = -p(0)/o(0)/o(x)Si(x,A,Ao). 
The SPPS summations Sq and Si have the following expressions 

oo oo 

Eo(x, A, Ao) = J]x(2")(x)(A - Ao)", Si(x, A, Ao) = Y,X'^'--'\x){X - Ao)"'^ , 

n=0 71=1 

where the coefficients X^'^^(x), X^^\x) are given by the recursive relations 

7;X('^-i)(0/o'(Oc?e foranoddn 
X(")(x) = <J (33) 
-/;x("-i)(0^^(^ foranevenn 
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foranoddn 

= <^ (34) 

[/;X("-i)(0/o'(0^^e foranevenn , 

which are identical to Eqs. (ED-dH]) unless for obvious sign changes. 

One can check by a straightforward calculation that the solutions /i and /2 fulfill 
the initial conditions f l27|l . Having obtained the fundamental system of solutions for 
any value of A, one can apply the construction fl28|l in order to obtain the Bloch 
solutions which become eigenfunctions for A being eigenvalues. 



4.4 Hill's discriminant in SPPS form 



We are ready now to write the Hill discriminant D{X) = /i(T, A) + f2{T, A) in a simple 
explicit form using the SPPS expressions of fi(T, A) and f2(T, A) in ( 132|) 

Joiuj Jo[l ) 

+ (/^(0)/o(T) - /o(0)/^(T))p(0)Si(T, A, Ao) . 

Finally, taking into account that fo{x) is a T-periodic function /o(0) = fo{T) and writ- 
ing the explicit expressions for So(^) A, Aq) and So(T, A, Aq) we obtain a representation 
for Hill's discriminant associated with (126 p 

oo 

D{X) = J2 (x^^'^KT) + X(2")(r)) (A - Ao)". (36) 



n=0 



Thus, only one particular nodeless and periodic solution /o(x, Ao) of ( 126|) is needed 
for constructing the associated Hill discriminant. We formulate the result ( 136|) as the 
following theorem: 

Theorem 5 Let Xq be the lowest eigenvalue of the periodic Sturm- Liouville problem 
^2D\j on the segment [0,T] with periodic boundary conditions and /o(a:;, Ao) be the cor- 
responding eigenf unction. Then the Hill discriminant for ^2B\j has the form where 
and are calculated according to / f^) and ([^]), and the series converges 

uniformly on any compact set of values of X. 

To illustrate the formula f l5B]) we consider a simple example. Let q{x) = 0, p{x) = 1 
in equation fl2^ . It is easy to see that the associated discriminant is D{X) = 2 cos -\/AT, 
from where we obtain Aq = and a corresponding non-trivial periodic solution is 
fo{x) = 1. Now making use of this solution we construct the discriminant by means of 
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the formula The coefficients X(2")(r) and X(2")(T) given by ([MD and take 
the form 

X(2")(r) = X(2")(T) = (-1)"^, n = 0,l, 2, ... . 
The substitution in (!36|) gives L>(A) = 2 cos Vat. 

4.5 Construction of the first eigenvalue Aq by the SPPS method 

Notice that in the expression f l35|) for -D(A) and in all reasonings previous to it we do 
not use the periodicity of the solution fo{x, Aq), therefore (|35|) and the whole procedure 
for obtaining it are valid for any A* such that there exists a corresponding solution 
/*(x,A*) which is bounded on [0,T] together with Such a solution A*) 

can be obtained in the following way 

A.,) = A,^) + if^,2{x, AJ (37) 

where /*,i(x, A*) and f*,2{x,K) are given by ( l30i) with A=k instead of Aq. For more 
details see [62]. The pair of the independent solutions /i(x, A) and /2(x, A) of (126!) 
given by (!32|) of course are independent of the choice of the solution fo{x, Aq), hence 
instead of /o(x,Ao) in ( l35l) one can take /*(a;, A*) given by ( !37l) . Thus, in terms of 
/*(x, A*) where A* is essentially arbitrary, -D(A) can be represented as a series in powers 
of (A - A.) 

^(^) = E (||y^^"^(^) + (38) 

+ (/:(o)/*(r) - /.(o)/:(r))p(o)x(2"+i)(T)^ (A - A.)^ 

Now the band edge Aq required for the formula (1361) can be calculated as a first zero 
of the expression -D(A) — 2 where -D(A) is given by (1381) . For the numerical purpose 
it can be useful to know the interval containing Aq. Since g is a bounded periodic 
function, there is a number A which satisfies the inequality q{x) > A Va; G R. It is 
known [33] that -D(A) > 2 for all A ^ A, therefore the lower estimate for Aq is the 
following 

Ao ^ min q{x). 

The upper bound can be obtained considering the Rayleigh quotient for periodic prob- 
lems |86] 

^ ^ lo {P jx) {x)f + q jx) {u{x)f) dx 
!o{'^{^)fdx 

where u{x) G [0,7"] is periodic with period T. The equality occurs if and only if 
u{x) is an eigenfunction corresponding to Aq. 
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4.6 Hill's discriminant of the supersymmetric (SUSY)-related 
equation 

In this subsection, we consider the SUSY partner equation of Eq. fl2^ and write down 
the SPPS form of its solutions. The latter allow us to prove the equality between the 
Hill discriminants of equation ( 126|1 and its SUSY-related Eq. ( HT]) . For various aspects 
of SUSY periodic problems, see [STJ [251 121] • 

The left-hand side of the equation (126!) can be factorized in the following way [87] 



(-4p^(x) + + fix), (39) 

where 4 means the x-derivative, the superpotential ^(x) is defined as follows <l>(x) = 
—p^ix)j^^rjl^. Using this factorization the coefficient g(x) can be expressed as 

q(x) = $^(x) - (^p^ix)^ix)^ + Xq. 
Introducing the following Darboux transformation 

{pHx)d, + $(x)) fix, X) = fix. A), (40) 



one obtains the equation supersymmetrically related to equation ( l26l) 

(pHx)d^ + <I)(x)) (^-d^p^ix) + $(x)) fix,X) = Xfix,X), 
which can be written as follows 

- 4(p(x)4/(x, A)) + qix)fix, X) = Xfix, X) , (41) 
where g(x) is the SUSY partner of the coefficient qix) given by 

qix) = qix) +2p^ix)^'ix) - p^(x)(p^(x))" . (42) 

It is worth noting that as ^(x) is a T-periodic function, the Darboux transformation 
assures the T-periodicity of qix). In addition, when p(x) is a constant, the SL coeffi- 
cient g is a quantum-mechanical potential, while g(x) is its Darboux counterpart also 
termed a supersymmetric partner in quantum mechanics. 

The pair of linearly independent solutions /i(x. A) and f2ix,X) of (HT]) can be 
obtained directly from the solutions ( !32l) by means of the Darboux transformation (l40ll . 
We additionally take the linear combinations in order that the solutions /i(x. A) and 
/2(x, A) satisfy the initial conditions /i(0. A) = /2(0, A) = 1 and /{(O, A) = /2(0, A) = 

f\ix, X) ^ 4^^So(x, A) + ^!Mrb^^^E,(x, A), (43) 
pH^)foix) (A-Ao)/o(0)p^(x)/o(x) 

/2(x, A) = Si(x, A). (44) 

(A-Ao)/o(0)pi(x)/o(x) 
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These two solutions allow us to write the expression for Hill's discriminant associ- 
ated to the equation (HTi) . that is -D(A) = /i(T, A) + /2(T, A). This requires the expres- 
sion of the derivative of f2{x, A) and evaluating it for x = T. In addition, one should 
notice that as the functions fo{x, Aq) and p{x) are T-periodic, i.e., /o(0, Aq) = foiT, Aq) 
andp(O) = p(T), then obviously, the functions /o(a;, Xo),p^{x) and [p2(a;)]' possess the 
same properties. The result is [55] 



D(A) = So(T,A) + So(T,A) 



[pi(x)]'U.=o-$(0) 



,(A-Ao)/o(0,Ao)p5(T)/o(T,Ao) 
[p^(x)]'U.=t/o(T, Ao) +p^(T)/^(T, Ao) 



(A-Ao)/o(0,Ao)p^(r)/2(T,Ao) 



Si(r,A). 



The substitution $(0) = —p^ (0) j°|o'^°| clearly shows that the expression in brackets 
vanishes leading to the simple formula 

oo 

5(A) = So(T, A) + So(r, A) = ^ (X(2")(r) + X(2")(T)) (A - Ao)" , 

n=0 

which is identical to (l36l) and therefore 



D{X)^D{X). (45) 

Thus, we can make the following statement: 

Theorem 6 Let Aq be the first eigenvalue of and fo{x, Aq) the corresponding T- 
periodic nodeless eigenf unction. Then the Darboux transformation with $(a;) = 
~P^{^) 'fl^(l'x°J) le.ads to a SUSY-related Eq. (^7]) with the preservation of the Hill dis- 
criminant, i.e., Eq. holds. 

From the identity of discriminants ( l45l) it is clear that Aq gives rise to a nodeless 
periodic solution /o(a;, Aq) of Eq. (14T]) . Taking A = Aq in ( H3|) and (jH]) we get this 
eigenfunction in the form /o(a;. A, 



p^(x)fo(xM) 

Notice that, the factorization method can be applied to Eq. fHT]) with the super- 
potential $i(x) = — p^(x) "|'|i^'^°j . In this case, we obtain the representation 



q = ^l{x)- [p^{x)^i{x)^ +Ao, 



which reduces to the equality ( H2l) if one notices the relationship $i(a;) = ip^ix))' — 
$(x). It can be also shown that q = q, where q = q{x) + 2p2 [x)^i{x) — p2[x){p^{x))" 
is the superpartner potential of q{x). Thus, the Darboux transformation fj40|) with the 
superpotential $i(x) applied to Eq. ( HTj) does not produce a different potential. 
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4.7 Numerical calculation of the eigenvalues based on Hill's 
discriminant in SPPS form 

As is well known, see e.g., [70], the zeros of the functions D{X)±2 represent eigenvalues 
of the corresponding operator. In this section, we show that besides other possible 
applications the representation (!36|) gives us an efficient tool for the calculation of the 
discrete spectrum of a periodic Sturm-Liouville operator. 

The first step of the numerical realization of the method consists in calculation 
of the minimal eigenvalue Aq by means of the procedure given in subsection 14.51 and 
subsequently in construction of the corresponding nodeless periodic solution /o(x, Aq) 
using formula f l3T]) . The next step of the algorithm is to compute the functions X^"^ 
and X*^"^ given by fl33|) and flMl) . respectively. This construction is based on the 
eigenfunction /o(x,Ao). Finally, by truncating the infinite series for -D(A) f l36|) we 
obtain a polynomial in A — Aq 

N 

Dn{X) = J2 (^^'"H^) + A:(2")(T)) (A - Ao)" (46) 

n=0 

N 

= 2 + ^ (x(2«)(T) + X(2")(r)) (A - Ao)". 

n=l 

The roots of the polynomials Di^{X) ±2 give us eigenvalues corresponding to Eq. fl2B]) 
with periodic and antiperiodic boundary conditions. 

As an example, we consider the Mathieu equation with the following coefficients 

p{x) = 1, q{x) = 2rcos2x . 

The algorithm was implemented in Matlab 2006. The recursive integration required for 
the construction of xi"\ X^"'\ X(") and X^") was done by representing the integrand 
through a cubic spline using the spapi routine with a division of the interval [0, T] into 
7000 subintervals and integrating using the fnint routine. Next, the zeros of DAr(A) ±2 
were calculated by means of the fnzeros routine. 

In Tables 4.1 and 4.2, the Mathieu eigenvalues were calculated employing the SPPS 
representation (|36|) for two values of the parameter r. For comparison the same eigen- 
values from the National Bureau of Standards (NBS) tables are also displayed [83] . 
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Table 4.1: A„ for the Mathieu Hamiltonian 




rf 1 

/ X 


rf 1 

/ X 


I V 




A fNRS "1 


n 

W 


—0 4^^1 390^^973837 


—0 4551 3860 


1 

X 


—0 1 1 0948490387377 


—0 1 1 094889 


2 


1.859107160521687 


1.85910807 


3 


3.917024962694820 


3.91702477 


4 


4.371299312651704 


4.37130098 


5 


9.047736927007582 


9.04773926 


6 


9.078369587941564 


9.07836885 


7 


16.033018848985410 


16.03297008 


8 


16.033785039658117 


16.03383234 


9 


25.020598536509114 


25.02084082 


10 


25.021087773318282 


25.02085434 




Table 4.2: A„ for the Mathieu Hamiltonian 




r = 5 


r = 5 


n 


A„ (SPPS) 


A„ (NBS ) 





-5.800045777242780 


-5.80004602 


1 


-5.790080596840196 


-5.79008060 


2 


1.858191484309548 


1.85818754 


3 


2.099460384254221 


2.09946045 


4 


7.449142541577460 


7.44910974 


5 


9.236327731534002 




6 


11.548906947651728 




7 


16.648219815375526 




8 


17.096668282587867 




9 


25.510753265631860 


25.51081605 


10 


25.551677357240167 


25.54997175 



Figures 1 and 2 display the plots of the calculated Hill discriminants for two values 
of the Mathieu parameter. A supersymmetric Mathieu potential can be written in 
terms of the even Mathieu cosine function as follows |74] : 

V. = 2 f 'h£f^) % 2A„ - 2r co=(2.) (47) 
V Ce(Ao,r,x) / 

and has the same Hill discriminants for identical values of the parameter r. 



As another example, consider the potential 



V^i = (1 - cos Ax) - 3^ cos 2x , (4J 
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Figure 1: The polynomial Di\f{X) for the Mathieu equation with the parameter r = 1 
calculated by means of formula for N = 100. 



Figure 2: Same as in the previous figure but for r = 5. The first minimum goes down 
to -292.0066. 
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which belongs to the quasi-exactly solvable family of the so-called trigonometric Razavy 
potentials [HH]- The parameter ^ is a positive real number. In Tables 4.3 and 4.4, the 
Razavy eigenvalues were calculated employing the SPPS representation f p6|) for two 
different values of the parameter ^. For comparison we use the eigenvalues given by 
Razavy analytically in terms of the parameter ^ as follows [89] 



Ao = 2 (l - x/l + e 



A, 



A, 



2(1 + v/r+e 



Table 4.3: A„ for the Razavy Hamiltonian 




e = l 


e = l 


n 


\n (SPPS) 


A„ (Ref. [89J ) 





-0.828430172936322 


-0.828427124746190 


1 


-0.627099642286704 




2 


2.315154289053194 




3 


3.999948252396118 


4 


4 


4.834668005757639 


4.828427124746190 


5 


9.246360795065604 




6 


9.305957668676312 




Table 4.4: A„ for the Razavy Hamiltonian 




e = 2 


e = 2 


n 


\n (SPPS) 


A„ (Ref. [89J) 





-2.472136690058546 


-2.472135954999580 


1 


-2.428288532265432 




2 


3.193559545313260 




3 


4.000042398350143 


4 


4 


6.472170127477180 


6.472135954999580 


5 


9.864070609921770 




6 


10.253303565368553 





In Figs. 3 and 4 we display the plots of the Hill discriminants for the values of 
the Razavy parameter ^ = 1 and ^ = 2, respectively. According to our results in 
subsection (14.61) . the Hill discriminant is the same in the case of the supersymmetric 
partner potential 



V2 = Vi+A cos 2x 



2^(0 



i-A{i) cos2x 



8A(0 sin^ 2x 
{i-A{i) cos2x)^ 



(49) 



for the same values of the parameter ^. In the latter equation A{^) = |^1 — a/1 + 

As final comments to this subsection, we believe that the calculation of the Hill 
discriminant through f H6|) offers clear numerical advantages with respect to other more 
complicated formulas for this important quantity provided in the literature, such as 
Jagerman's so-called cardinal series representation [50], the infinite determinant rep- 
resentation involving the Fourier coefficients of the potential as well as the spectral 
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Figure 3: The polynomial Dn{X) for the Razavy equation with the paranieter ^ = 1 
calculated by means of formula for N = 100. 

parameter in the book of Magnus and Winkler [76], a matrix representation whose 
entries are comphcated phase integrals obtained by the phase-integral method used 
by Froman [30], and Boumenir's representation in terms of integrals derived from the 
inverse spectral theory [H]. 

5 Spectral and transmission problems on the whole 
line 

In this section we consider the one-dimensional Schrodinger equation 

Hu{x) = —u"{x) + Q{x)u{x) = \u{x), x G M, (50) 

where 

{«!, X < 0, 

q{x), 0<x<h, (51) 
a2, X > h, 

ai and 02 are complex constants and g is a continuous complex-valued function defined 
on the segment [0, h]. Thus, outside a finite segment the potential Q admits constant 
values, and at the end points of the segment the potential may have discontinuities. 
We are interested in two classical problems. The first is the quantum-mechanical 
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Figure 4: Same as in the previous Ggure but for ^ = 2. The Grst minimum of Hill's 
discriminant goes down to -55.01. 

spectral problem, we are looking for such values of the spectral parameter A G C for 
which the Schrodinger equation possesses a solution u belonging to the Sobolev space 
H'^{R) which in the case of the potential of the form (IHTj) means that we are looking 
for solutions exponentially decreasing at ±oo. 

The second consists in finding the reflectance and transmittance of the inhomoge- 
neous layer described by q. We will formulate this problem in the form in which it 
arises in electromagnetic theory though both problems come not from one but from 
many different branches of physics and engineering. 

5.1 Quantum-mechanical spectral problem 

The eigenvalue problem f l50|) is one of the central in quantum mechanics for which 
H is a self-adjoint operator in L^(M) with the domain if^(M). It implies that Q 
is a real-valued function. In this case the operator H has a continuous spectrum 
[min {«!, 0:2} , +00) and a discrete spectrum located on the set 

[ min g(a;), min {«!, 0:2}). (52) 

xe[o,h] 

Computation of energy levels of a quantum well described by the potential Q is a 
problem of physics of semiconductor nanostructures (see, e.g., Other important 

models which reduce to the spectral problem ( l50l) arise in studying the electromagnetic 
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and acoustic wave propagation in inhomogeneous waveguides (see for instance [1], [2T], 

m, m, m, m, m)- 

Hence in the applied problems it is important to have effective and rapid numerical 
methods for the solution of the problem (1501) . The most frequently applied is the 
shooting method (see, e.g., It has well known limitations due to the intrinsic 

difficulties of the shooting procedure, especially when the spectral parameter as in the 
problem under consideration participates in the boundary conditions (see equalities 
( 154|) and (1571) below). It is much more convenient to have available an analytical form 
of a dispersion equation associated with the eigenvalue problem. In that case solutions 
of the dispersion equation can be approximated using different numerical techniques. 
However the dispersion equation is available only in really few examples (see [55]). 
There is another method developed for symmetric potentials |13]. Below we compare 
numerical results of our approach with the results from ^3]. 

For simplicity we will assume that ai and a2 are real constants and the function 
g is a continuous real- valued function on [0, h] though the presented method is appli- 
cable to the more general situation when Q is complex valued. In this case necessary 
modifications must be made mainly in the reduction of the original problem on the 
whole line to a problem for the equation 

— u" (x) + q{x)u{x) = Xu{x) , xE{0,h) (53) 

with three boundary conditions at the end points of the interval (0, h) (see below) 
meanwhile the application of the SPPS method suffers no essential changes. Our 
analysis follows that from [T8] . 

For X < we have to consider the equation —u" + («i — A)^ = 0. Its solutions 
decreasing at — oo exist if only ai — A > 0. Denote /i = +^/a^'^^. Then the required 
solution for x < has the form u{x) = cic^^ and the multiplicative constant can 
always be chosen equal to one. Thus, u{x) = e^^, a; < 0, from where 

m(0) = 1 and u{0) = /i. (54) 

This gives us the initial conditions for the solution on the interval (0, h) which we 
will construct following Theorem [H For that we need first a nonvanishing particular 
solution of the equation 

- u'o'(x) + q{x)uo{x) = (55) 

which as was explained in Section [2] can be constructed by means of the same SPPS 
method. Indeed, formulas f|T4l) - f|T7j) where p should be chosen equal to —1 and xq = 
give us a couple of linearly independent real-valued particular solutions vi and V2 of 
fl55|) . Hence (see Remark H]) the required nonvanishing solution of fl55l) can be chosen 
as uq = vi + iv2- Let us notice that as mo(0) = 1 and Uq{0) = i the initial conditions 
satisfied by the solutions of ( l53l) ui and U2 constructed according to ([5]) have the form 

Mi(0) = Mo(0) = 1, u[{0) =u'o{0) =t, 
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n2(0) = 0, n'2(0) = -— - = -1. 

From these relations we obtain that the solution of f lS^ satisfying the initial conditions 
( 154|) has the form 

u{x) = ui{x) + {i — fi)u2{x) < X < h. (56) 
In the region x > h the solution of equation fjSOj) has the form 



u 



from which we obtain that the existence of an eigenfunction is possible if only 1/0:2 — A G 
M. Hence 0:2 > A and we denote v = +1/02 — A. Consequently, = Ce~'^^^~^^ and 
= C, u'{h) = —vC where C is an arbitrary constant. Thus, the eigenvalues of 
the problem are such values of A for which the solution ( l56l) satisfies the condition 

u\h) + uu{h) = (57) 



where, as above, u = +1/0:2 — A and 02 > A. 

In order to write down the explicit form of the dispersion equation ( 157|) in terms 
of the spectral parameter power series we calculate the derivatives of the solutions of 



/ \ o^ / -| o^ 

^; = _ Ay A"X(2"+i) and = -W2 - -yA"X(2-). 
Thus the derivative of the solution f l56|l has the form 

u' = ^u--l ^A"+ix(2"+i) + (z - /x)5^A"X(2") 

\n=0 n=0 / 

Substituting this expression into ( IFT]) we arrive at the following result obtained in 
and formulated here in the form of a theorem. 

Theorem 7 Letai, 02 &e real numbers, q be a real-valued continuous function defined 
on [0,h] and Q be defined by 137]) . Then A G [mina.g[o,h] Q'(a;), min {oi, 02}) is an 
eigenvalue of the problem / fJOj) z/ anc? on/y if the following dispersion equation 



(00 00 

^A"X(2")(/i) + - v^^T^) ^A"X(2"+i)(/i) 
ra=0 n=0 
_ / 00 00 

|- 5^A"+ix(2"+i)(/i) + (z - v^^T^) ^A"X(2")(/,) 

' \n=0 n=0 



/ 00 00 \ 

+ v/^^no(/i) 5^A"X(2")(/i) + (^z - ^A"X(2"+i)(/,) = , (5^ 



vn=0 ra=0 
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is satisfied and the corresponding (unique up to a multiplicative constant) eigenf unction 
has the form 

e^^, a; < 0, 

u{x) = { ui{x) + {i — fi)u2{x), < X < h, 

(mi(/i) + (i-/i)M2(/i))e-''(^-'*), x>h, 



where n = +y/ai — X, v = +\^a2 — A and ui, U2 are defined by ([^ where uo is the 
nonvanishing solution of 153]) on {0,h) satisfying the initial conditions Uo{0) = 1 and 
Mq(0) = i, p = —1, r = 1 and Xq = 0. 

All the coefficients in equation ( l58l) : uo{h), Uq^H), X^'^\h) and X^'^\h) are easily 
and (as our numerical tests show) accurately obtained from the definitions introduced 
above, and the roots of the dispersion equation coincide with the eigenvalues of the 
problem and can be found using many available methods. 

In what follows, let us consider a relatively simple situation: ai = a2- Rearranging 
the terms in equation (150|1 this case always can be reduced to the case cti = ^2 = 0. 
Then u = fx = \/—\, jj? = —A and A" = (— The dispersion equation takes the 
form (here we correct some easily detectable misprints in [18]) 

u',{h){l + iX^'\h)) 



Uo{h) 

+5^(-i)V"K(/i)^^'"H/^) — --x^^''-'\h) + iu'oih)x'^^''+'\h) 

n=l 

^X(2")(/i)+«o(/i)X(2"-i)(/,)) 

Uo{h) 

oo 

+ ^(-l)V™+'(-Wo(/^)^^'"+'n^) + — TTT^^""^/^) 
m=0 > 

+«Uo(/i)X(2-+i)(/i) + uo(/i)X(2™)(/i)) = 0. 
Thus, the dispersion equation has the form 

oo 

J2 = (59) 

fc=0 

where 

a, = u',ih) {1 + tX^'\h)) - ^ (60) 

a2n = {-inu',{h)X('-\h) - ^x(2-i)(/^) + z«;,(/,)x(2"+i)(/,) 

uo{h) 

' X(2")(/i) + uo(/i)X(2"-i)(/i)), riGN, (61) 



uo{h) 
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a2n+i = (-l)"(-w'o(/i)X(2"+i)(/i) + _-X(2-)(/i) 

Uo{h) 

+mo(/i)X(2'»+i)(/i) n = 0,1,2,.... (62) 

The problem is reduced to the problem of finding zeros of an analytic function given 
by its Taylor series with the coefficients a^, k = 0,1,2, .. .. 

The usual approach to numerical solution of the considered eigenvalue problem 
consists in applying the shooting method (see, e.g., [S]) which is known to be unsta- 
ble, relatively slow and to the difference of our approach does not offer any explicit 
equation for determining eigenvalues and eigenfunctions. In |13] another method based 
on approximation of the potential by square wells was proposed. It is limited to the 
case of symmetric potentials. The approach based on the SPPS is completely different 
and does not require any shooting procedure, approximation of the potential or nu- 
merical differentiation. Derived from the exact dispersion equation (159|) we consider its 
approximation X]^o '^fc/^*^ = and in fact look for zeros of the polynomial XlfcLo '^kf^'^ 
in the interval [ming(x), 0). Here we give only one example of numerical computation 
of eigenvalues referring to [18] for more examples and discussion. 

We consider the potential Q defined by the expression Q{x) = — t;sech^x, x G 
(— oo, oo). It is not of a finite support, nevertheless its absolute value decreases rapidly 
when X — )■ ±oo. We approximate the original problem by a problem with a finite 
support potential Q defined by the equality 

X < —a 
a < X < a 
X > a . 

An attractive feature of the potential Q is that its eigenvalues can be calculated 
explicitly (see, e.g., [39]). In particular, for v = m(m + 1) the eigenvalue A„ is given 
by the formula A„ = — (m — n)"^, n = 0,1, . . .. 

The results of application of the SPPS method for v = 12 are given in Table 15.11 
in comparison with the exact values and the results from [43j. 



Table 5.1: Approximations of A„ of the Hamiltonian H = —D^ — 12sech^x 


n 


Exact values 


Numerical results from [431 


Numerical results using SPPS {N = 180) 





-9 


-9.094 


-8.999628656 


1 


-4 


-4.295 


-3.999998053 


2 


-1 


-0.885 


-0.999927816 



Qix) 



0, 

-V sech^ 
0, 
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5.2 Transmission problem for inhomogeneous layers 

In this subsection we apply the SPPS method to the problem of finding the refiectance 
and transmittance of a finite inhomogeneous layer. This is a classical problem which 
still attracts a lot of attention due to its numerous applications in modern engineering, 
optical physics, solution of nonlinear problems and many other fields. Different meth- 
ods for numerical solution of the problem have been proposed, mainly based on well 
known canonical techniques for approximate solution of ordinary differential equations 
such as the finite differences or expansion in power series (see, e.g., [56], |101] . [5^ . 
[19]). One of the most used methods involves the approximation of the inhomogeneous 
layer by a structure consisting of many homogeneous layers (see, e.g., [50], [SI], [2Z], 
[96j). Asymptotic methods such as the perturbation method or the WKB method are 
also applied to this problem (see, e.g., [12], [SZ], [H2], |101] )- though in the case of a 
finite inhomogeneous layer the WKB technique does not seem advantageous. Mean- 
while the mentioned numerical approaches can give satisfactory results for certain 
fixed parameters of the problem their applicability is questionable when the solution 
of the problem is required, for example, for many different angles of incidence. The 
treatment of the oblique incidence case is not only interesting because of the many ap- 
plications in which that incidence is needed - in optical filters, light couplers -, but also 
because sometimes the interfaces are rough - their effects and analysis depending on 
their size -, and/or are not parallel (see, e.g., [82]). This is due to imperfect deposition 
conditions. Such problems in the generation of the inhomogeneous layer (or multi- 
layer) have generated systems in which the feedback of a reflectance, transmittance, 
or scattered light measurement is used to characterize the layer as it is created and to 
correct any discrepancies with the pre-established values. Such application must be 
able to recalculate the required correction proflle and requires a real-time computation 
of transmittance and reflectance. 

The mathematical statement of the problem involves a Helmholtz equation with a 
coefficient which is an arbitrary continuous function on a flnite segment and constant 
outside. More precisely, the scalar function u which represents a component of a 
linearly polarized electromagnetic wave in the case of an s-polarization satisfles the 
Helmholtz equation 



where u stands for the transverse component of the electric fleld, and in the case of a 
p-polarization satisfles the following Sturm-Liouville equation (see, e.g., [19], [35] ) 



in which v represents the transverse component of the magnetic fleld. Here k is the 
free-space circular wave number. The refractive index n preserves constant values ni 
and n2 in the regions x < and x > d respectively and is an arbitrary continuous 
function in the interval < x < d (see flgure [5]). For simplicity we assume n to be 




(63) 




(64) 
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Figure 5: An inhomogeneous layer. 



real valued though the method is equally applicable to the case of a complex refractive 
index. 

The propagation constant /3 is related to the angle of incidence of the wave in 
the following way /3 = fcsin6' (see, e.g., [19]), and /3 vanishes in the case of normal 
incidence. 

In spite of the fact that equations (j63|) and flMj) describe the behaviour of different 
components of an electromagnetic wave, corresponding to an electric and a magnetic 
field respectively, there exists a simple transformation from flM|l to fl63|l and vice versa 
(see, e.g., ^49j). Namely, if w is a solution of flM|) then U = v/n is a solution of the 
equation 

U"ix) + [eN^ix)~P^]Uix) =0 

where k'^N'^ = k'^n^ + n" /n — 2 [n' /nf' . Thus, in both cases the problem reduces to 
an equation of the form ( |63l) . 

We denote ki = -^Z k'^nf — P'^ and k2 = \/k^r^ — The solution u of fl63|l or v 
of flM|) respectively together with their first derivatives must be continuous at all x 
including the points x = and x = d. The incident wave in the region I (see figure 
E]) is assumed to have the form e~*'^i^, and together with the reflected wave the whole 
solution for a; < is the combination 

u{x) = e-'^''' + i^e'^^i", X < 

where the constant R is the reflection coefficient whose absolute value is less than 1. 
The solution corresponding to the transmitted wave in the region II has the form 

u{x) = Te''^^'', x>d 

where T is the transmission coefficient. In the case of unabsorbent media for the 
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Figure 6: Incident, reflected and transmitted waves. 



normally incident waves the following energy conservation relation holds 

\Rf + n2 \Tf /rii = 1. 



(65) 



Let us suppose that the two linearly independent solutions yi and 1/2 of in 
the interval of inhomogenicity < x < d are known such that the following initial 
conditions are satisfied: 

yi(0) = l, y[{0) = Q (66) 

and 

1/2(0) = 0, 1/^(0) = 1. (67) 

Then we are able to obtain analytic expressions for R and T in terms of Ui and U2 
[H]. We have 

-kik2y2{d) - y'i{d) - ik2yi{d) + ikiy'^id) 



and 



R 



T 



[y'i{d) - kik2y2{d)] + i[k2yi{d) + kiy'^id )] 
2iki[yi{d)y'2{d) - y[{d)y2\ 



(6J 



(69) 



[y[{d) - kik2y2{d)] + i[k2yi{d) + kiy'^id)] ' 

These formulas remain valid for equation when one substitutes yi and y2 with the 
solutions vi and V2 of flMl) satisfying the initial conditions (!66|) and (!67|) respectively). 

Thus, the transmission problem for an inhomogeneous layer consists in computing 
a couple of solutions of fl63p (or (1641) ) in the interval of inhomogenicity < x < (i, 
satisfying the initial conditions (16 6 p and (16 7p . and then the reflection and transmission 
coefficients are found from fl68l) and fl69l) . For computation of these solutions we use 
theorem [1] and take into account ffTTl) and f|T2|) where it is convenient to choose xq = 0. 
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There are several examples of explicitly solvable inhomogeneous profiles (HD] , |101] • 
These were used in [T7] for testing the results obtained by means of SPPS. In all 
numerical simulations the achieved accuracy was remarkable. 

6 Zakharov-Shabat eigenvalue problem 

In this section we study the Zakharov-Shabat system with a real-valued potential. It 
arises in the solution via the inverse scattering method of several nonlinear evolution 
equations such as the nonlinear Schrodinger equation, the sine-Gordon equation and 
the modified Korteweg-de Vries equation. For example, in the case of the nonlinear 
Schrodinger equation, eigenvalues of the Zakharov-Shabat system correspond to soliton 
solutions implemented in fiber optics (see, e.g., [IS]). The assumption that the poten- 
tial is real valued is natural and common in the engineering literature,- it includes the 
conventional profiles such as the rectangular, the Gaussian and the hyperbolic secant. 

In [68] a general solution of the Zakharov-Shabat system with a real potential in 
terms of SPPS was obtained and used for deriving a dispersion equation corresponding 
to the eigenvalue problem with a compactly supported potential. Once again the 
problem is reduced to a problem of localizing zeros of an analytic function given by its 
Taylor series. For numerical approximation of eigenvalues one can consider a truncated 
series and thus for practical computation the eigenvalue problem reduces to finding 
roots of a polynomial. 

The Zakharov-Shabat system with a real potential has the form jlU2j . [70] 

dni{x) — Xni{x) = U{x)n2{x) (70) 
dn2{x) + \n2{x) = —U{x)ni{x) , (71) 

where 9 := f/ : M — ?■ M is the potential and U E Li(— oo, oo); the solutions rii and 
n2 in general are complex valued and the spectral parameter A is a complex constant. 
It is convenient to rewrite the Zakharov-Shabat system using the following notations 

u = 111 + in2, V = rii — in2, q = iU . 
Then (ITOjl . (ITTjl takes the form of a Dirac system with a scalar potential (see, e.g., [T5] . 

HZ], [H], M) 

{d + q{x))u = Xv, (72) 

{d — q{x))v = Xu . (73) 

From these equalities it is easy to see that u and v are solutions of the following 
second-order differential equations 

{d - q{x)){d + q{x))u{x) = X\{x) (74) 

and 
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{d + q{x)){d — q{x))v{x) = A^f (x) . (75) 

The differential operators on the left-hand side can be written in the form of stationary 
Schrodinger operators describing supersymmetric partners 

{d - q){d + q) = d'^ + {dq - q^) and {d + q){d - q) = d'^ - {dq + q"^). 

Nevertheless, precisely the factorized form ([71]), fl75]) presents certain advantage for 
applying the SPPS method due to the possibility to write down closed-form solutions 
of dHD and ([75]) for A = 0. Namely, let Q{x) = J q{x)dx. Then uo{x) = e'^^^) and 
vq{x) = e'^^^^ are solutions of the equations {d—q){d+q)uo = and {d+q){d—q)vo = 
respectively. Note that for a continuous function q defined on a closed finite interval 
both uq and Vq are devoid of zeros. 

The systems of auxiliary functions {X^"^^^_^ and < X*^") > in this case are de- 
fined as follows 

X(°)(x) =X(°)(a;) = 1, (76) 



X(")(x)=/ X("-i)(s)e(-^)"2Q(.)^^^ (77) 

Jxo 

X 

X(")(x) = / X("-i)(s)e(-i)"^'2Q(s)^^_ (7g) 



We obtain the following SPPS form of a general solution of the Zakharov-Shabat 
system. 

Theorem 8 fUE/ Let U be a continuous real-valued function defined on a finite segment 
[a, b] C M. Then the general solution of the Zakharov-Shabat system [7U\ ), [71\ ) has the 
form 

oo oo 

n,(x) = £1 Ve(-i)"«(^')A"X(")(x) + ^ V e(-i)"^'«(^)A"X(")(x), (79) 
2 ^ 2A ^ 

n=0 n=0 
oo . oo 

^2(x) = ^ V(-I)'^e(-1)"«(^)A"X(")(X) - ^ V(-l)"e(-^)"^'<3(x);,n^{n)(^)^ (gg) 

2 ^ 2A ^ 

n=0 n=0 

where Ci and C2 are arbitrary complex constants, Q is an antiderivative of q = ill, 
Xq G [a,b] and the series converge uniformly in [a,b]. 

Solutions of the Zakharov-Shabat system (!70l) . (!7T]) satisfying the following asymp- 
totic relations 

a(x,A) ^ Q^e^"^, x-^-cx), (81) 
e(x,A)- (°)e-'^ x-^oo (82) 
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are called the Jost solutions [T]. The eigenvalue problem for the Zakharov-Shabat 
system with a real valued potential consists in finding such values of the spectral 
parameter A for which Re A > and there exists a nontrivial solution n satisfying the 
Jost conditions (EI]) and (I82]). 

If the real valued potential U has a compact support on the segment [—a, a] it is 
easy to see that the eigenvalue problem reduces to find such values of A (Re A > 0) for 
which there exists a solution of flTO]) . f lTT]) satisfying the following boundary conditions 



ni{-a) 
ni{a) 



1, 
0. 



n2{-a) = 0, 



^3) 



We refer here to [58] and [59] for estimates of the number of real eigenvalues of a 
compactly supported potential. 

The next statement gives us a dispersion equation equivalent to the Zakharov- 
Shabat eigenvalue problem for the real, compactly supported potentials. 



Theorem 9 JUB^ Let U be a continuous real-valued function with a compact support 
on the segment [—a, a]. Then A (He A > 0) is an eigenvalue of the Zakharov-Shabat 
system if and only if the following equation is satisfied 



oo 

^ A" (e(-i)"«('^)x(")(a) + e(-i)"^'«('^)x(")(a)) = 



^5) 



n=0 



where Q{x) = i J U{t)dt and Xq = —a in ^-(7^. 

—a 

If A is an eigenvalue then the corresponding eigenvector is given by 

ip + (p 



It 



with 



ip{x) 



if{x) 



i^iix] 

1p2{x] 



/ 1 oo 



n=0 



1 ^(-l)«e(-i)"«(^)A"X(")(x) 

\ 2 n=0 



/ 1 oo \ 

( A ^ e(-i)"^'«(^)A"X(")(x) \ 

2 n=0 



v 



^(_l)™e(-ir+^Q(-)A"X(")(x) 

2 n=0 



17) 



The theorem reduces the Zakharov-Shabat eigenvalue problem with a compactly 
supported potential to the problem of localizing zeros (in the right half-plane) of an 

oo 

analytic function k(A) = ^ cinA" of the complex variable A with the Taylor coefficients 

n=0 

ttn given by the expressions 
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Equation (15^ represents a dispersion equation of the eigenvalue problem. The coeffi- 
cients be easily and accurately calculated following the definitions introduced 
above. For the numerical solution of the eigenvalue problem one can consider a poly- 
nomial 

N 

'^7v(A) = ^a„A" (89) 

n=0 

approximating the function k. For a reasonably large N its roots give an accurate 
approximation of the eigenvalues of the problem. 

As a numerical example we consider the rectangular box referring the reader to [SB] 
for further numerical tests and details. For the rectangular box the exact solution sat- 
isfying the boundary conditions (!83|) is known. Such system can be applied to describe 
the problem of the diffraction of a wave by a screen with a slit (see [75]). Discrete 
eigenvalues of the spectral parameter can also be approximated using a variational 
principle approach [S2]. The potential is defined by the equality 

U^^^ = / ^' 1^1 ^ 

\ 0, elsewhere. 

A dispersion equation in this case can be obtained explicitly and written as follows 

7 cos 27 + A sin 27 = , (90) 



where 7 = v^A^ — A^. 

For solving the dispersion equation (!90|) the routine NSolve of Wolfram Mathe- 
matica 7 was used. We considered a = 1. In the case A = 1 the routine NSolve 
delivers one solution of (ED]) A^^°'^^ = 0.31902252414261895. Apphcation of the SPPS 
method with m = 2000 and = 120 gives us the value Aq = 0.31902252414254. The 
agreement is up to the 12th digit. Taking m = 4000 and A^ = 180 we obtain still a 
better approximation Aq = 0.319022524142619. The agreement is up to the 14th digit. 

In the case A = 4 there are three eigenvalues. NSolve delivers the following val- 
ues A^^"''''^ = 0.41262411401896715, X^^oive ^ 2.8945478628320327 and A^^'''^'^ = 
3.749624961605374. Application of the SPPS method with m = 2000 and A^ = 
100 gives us the values Aq = 0.412624114002, Ai = 2.8945478628329 and A2 = 
3.7496249616095, and with m = 4000 and A^ = 180: Aq = 0.4126241140179, Ai = 
2.89454786283226 and A2 = 3.7496249616045. 



7 Conclusions 

We presented a review of recent research and applications of spectral parameter powers 
series (SPPS) representations for solving initial and boundary value problems as well as 
spectral and related problems for Sturm-Liouville equations. Application of the SPPS 
approach allows one to obtain explicit analytic forms of characteristic equations for a 
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variety of problems. Approximation of these equations represents a powerful, universal 

and accurate numerical method highly competitive with the best purely computational 
techniques. The SPPS method is algorithmically simple and can be easily implemented 
using available routines of such environments for scientific computing as Matlab. 
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